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ABSTRACT 


The inference of atmospheric structure 
from satellite radiometric observations re- 
quires an inversion algorithm. The ideal 
inversion technique should be accurate, self- 
limiting, free from bias, stable against 
noise, felxible, and simple in application. 

A variety of techniques has been spawned to 
meet these demands. One class, the nonlinear 
inversion methods, copes with the serious 
problem of data noise in an unusual fashion. 
Unlike linear techniques which require a 
priori data smoothing, the nonlinear method 
can be applied directly to raw data. The 
algorithm discriminates the noise input by 
resolving the inferences into two types of 
solution, associating the real roots with 
atmospheric structure while ascribing the 
Imaginary roots to noise. 

The algorithm appears capable of fur- 
ther refinement with the possibility of in- 
ferring systematic noise structure, i.e. 
pinpointing a channel consistently in error. 
An example is given of this error-sensing 
ability of the nonlinear inversion technique. 


* On leave from GCA Corporation, Bedford, Massachusetts 01730 
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1. INTRODUCTION 

The advent. of meteorological satellites has led to intensive investigation 
of the possibilities and limitations inherent in remote atmospheric sensing. 

The upwelling radiation intercepted by the orbiting spectrometer is an integral 
transform of the physical state of the atmosphere, symbolically expressible as 

I[*( * )] = Q { B[T(u)]}, (1) 

with x the monochromatic absorption coefficient at frequency v and B the 
Planck intensity, an implicit function of the vertical absorber distribution 
u through the temperature. For nadir frequency scanning in the far infrared 
the integral operator follows from radiative transfer theory as a Laplace* 
transform 

oo 

1(0.) l/x) = jf B(u)e xdu = xL[B(u)]. (2) 

Clearly the deduction of vertical thermal structure requires the solution of 
the inverse problem 



. might appear at . first sight as if we were making a problem where none 

existed. Certainly the inverse Laplace transform is known for a variety of 
functions and extensive tabulations are readily accessible. Two features 
both related to the observation of real data, intrude into the pure functional 
space of the mathematician. These are the unavoidable noise residing in the 
data and the finite. number of sensing channels. The fact of noise is partic- 
ularly, pernicious, imposing a sharp upper limit to the information deduction. 
Inversion, akin to differentiation, tends to amplify data error. The finite 
channel number forces one to infer a continuous vertical structure from a 
united sampling of the radiation field. Thus the goal of inversion cannot 
be to infer all the temperature structure. This is inaccessible in principle 
because of the presence of noise. Our aim must be the more modest one of 
seeking the optimum inversion technique which yields all the valid inferences 
inherent in the. observations. A more subtle corollary to this is that the 
inversion algorithm must be capable of discriminating noise from the data, 
at is a proper attribution of the signals either to radiating atmospheric 
sources or to extraneous noise. 


Let us set forth the requirements of an optimum, ideal inversion proce' 
dure c 

l) Accuracy: The inversion should reconstruct all the temperature 

structure inherent in the data. 


,, x? Self ' -limiting : The inversion method should not infer more structure 

than that permitted by the noise level. 

, 3) Freedom from bias: The inference should not be weighted towards any 

predetermined structure, such as a climatological set. Neither should the 
functional representation chosen force an unnatural configuration on the 
thermal profile. 
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4) Noise stability: The algorithm .should be capable of discriminating 

noise from signals arising from emitting atmospheric sources. 

5) Flexibility: The technique should be applicable to any meteorological 

situation and to any choice of channel sensing input. 

6) Simplicity: The inversion procedure should be compatible for com- 

puter programming and real time temperature data readout. 

It is a pleasure to report that a nonlinear inversion method has been 
developed which gives promise not only of fulfilling these stringent demands, 
but appears capable of discriminating between systematic and random errors. 

That is to say the inversion algorithm will go beyond indicating the presence 
of noise, by pinpointing which channel is in error and by how much! To indi- 
cate how this can be done, let us review the history of the inversion problem 
in sufficient depth to place this nonlinear Fourier inversion method in per- 
spective. 


2. INVERSION THEORY AND CRITIQUE 

Most inversion methods are basically linear techniques modified in 
various ways to cope with the noise problem. They proceed by expanding the 
Planck intensity in a suitably chosen orthogonal polynomial set 

n 

B(u) = 2 B.P .(u). (4) 

j=l J J 

Substituting this expression into Eq. (2) leads to a linear simultaneous set 
of equations which must be satisfied by the n measured intensities 

n 

I(0,1/k^) — 2 Bjpj(xq) + e(x^) i 1, 2,,..,n, (5) 

j=l 

where p-j(x) =xL[P.(u) and e is the noise vector. Experience has shown that 

the best choice for representing the source is a set of empirical orthogonal 
functions based on climatological data. The solution consists of a matrix 
inversion to determine the n weights of the prechosen orthogonal members, 
with the data appropriately smoothed by the subtraction of the error vector e. 

Nonlinear inversion is a completely different approach to the problem. 
Rather than specifying in advance the functional representation, the weights 
and members of the set are inferred directly from the intensity data. The 
data points are used to generate a unique characteristic equation whose 
eigenfunction solutions form the members of the set. 

In the first application to nonlinear inversion, the Planck intensity 
was approximated by spline functions, i.e., a series of slabs (step functions) 
or ramps. Since the nonlinear inversion method using spline functions is 
documented elsewhere (King 1964), we shall merely epitomize its merits and 
demerits relative to linear methods before proceeding with the new approach. 
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Three advantages come to mind. First is the fact that the intensity data 
determine directly the choice of members of the set. For spline functions 
this corresponds to the slab thicknesses or the distance between successive 
ramps. The second feature, the uniqueness of the inferred profile, follows 
as a corollary of the first. The thicknesses are given by the algorithm as 
the n roots, necessarily unique, of an nth degree polynomial. Stated ex- 
plicitly, for any suitably chosen 2n intensities, there is one and only one 
array of n slabs which will fit the data. The third feature is the most 
subtle, most unexpected, and most important. This is the response of the 
algorithm to noise in the data. In this event one or more of the roots be- 
come negative. These inadmissible roots characteristically have small weights 
associated with them. The remaining valid roots are relatively unperturbed 
and preserve a high fidelity representation of the temperature profile. Thus 
the algorithm acts as a filter, discriminating between the rapidly varying 
noise components and the lower frequencies associated with temperature 
structure . 

Balancing somewhat these three advantages of member choice, uniqueness, 
and noise discrimination vis-a-vis the linear methods are four restrictions 
associated with the nonlinear spline function inversion. First is the stipu- 
lation that the Planck intensity be approximated by a slab or ramp solution. 
For example, a constant tropospheric lapse-rate cannot be satisfactorily fit 
by a single ramp configuration. More serious, perhaps is the algorithm re- 
quirement that the channels be chosen at consecutive integral multiples of 
the absorption coefficient of the most transparent channel. The channel 
positions in practice are chosen out of engineering considerations, and it is 
highly unlikely that a choice on that basis would be optimal for the algorithm. 
A third condition is the need for an independent determination of the tempera- 
ture at the top of the atmosphere [B(0)j in the slab algorithm. For the ramp 
solution, inputs for both B(0) and B'(0) are necessary. Finally, the Prony 
algorithm is applicable only for transmitt ances of an exponential function 
type. 


Certainly the most interesting and potentially the most useful feature 
of the nonlinear approach lies in its treatment of noise in the data. A 
natural question is whether the nonlinear technique has an underlying physical 
basis in transfer theory or if the inadmissible roots are mere artifacts of a 
solution algorithm. Accordingly, an effort was made to establish the con- 
straints necessary to yield slab or ramp solutions to the transfer equation. 
The search proved fruitless. Indeed it appeared that solutions of the spline 
function type were incompatible with classical transfer theory. 

The pursuit of this logic has led to the formulation of a new wave theory 
of radiative transfer which contrasts to the corpuscular approach of classical 
transfer theory. Although the preliminary outlines are clear, at this date 
(July 1968) the theory is not yet complete. The paper, "Remote Sensing and 
Inversion Techniques: State of Art, Kine ( 1967 ), indicates the progress and 

general direction of this research. 

Our prime concern here is not transfer theory itself, but rather the 
insight it provides for inversion. In this vein it is a pleasure to report 
that the wave transfer theory has as a direct consequence a new inversion 
method based on nonlinear Fourier analysis. Let us back up a bit to see the 
problem and how the new inversion technique fulfills the need. 
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In the slab formulation the upwelling intensity is approximated by a set 
of exponentially weighted step functions of height Bj and thickness Uj+q -Uj. 


oo -xu 

l(0,l/x) - B(0) - S q e dB(u) 


-XUj 


= Z B ± e = 2 Bpc,* , 

J J 


(6) 


where we have substituted u- = -in x^. By specifying integral values of 

J 

we are led to a nonlinear simultaneous equation set which is mathematically 
equivalent to the moment problem in physics. This set 


1(0, lA ± ) - B(o) = a ± 


B j X j 


— 0 , 1 , . . . , 2n— 1 


(7) 


possesses a unique solution of n slabs which can be obtained using the Prony 
algorithm. 

We are led to ask the three questions 

1) Are there other nonlinear sets soluble by a similar algorithm? 

2) Are these sets more flexible, i.e., less rigid than the slab or ramp 
configurations? 

3) Can the intensity sampling points be arbitrarily chosen, free from 
the Prony requirement of equally spaced intervals? 

The answer to these questions is affirmative, leading to the hope that we 
are approaching the goal described earlier of an optimum inversion method. 


3. NONLINEAR FOURIER INVERSION 

3.1 Application to Noise-Free Data 

We begin by noting that the most general solution of the transfer equation 
involves waves, i.e., exponential functions of imaginary argument, as well as 
the attenuating exponentials of classical theory 


I(u,l/x) = c 



b .e 
3 


-iw.u 

J 


+ 


x + i-us 


+ _1 + u + Q 
x 


(a) 
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By specifying odd parity for the Planck intensity B(u), the upwelling 
intensity l( 0 ,lA) and its inverse Laplace transform B(u) become 


1(0, iA) = c 


oo 

V 


J=1 


* b w 
j 3 


2 2 
x + Wj 







B(u) = L 1 

1 

(H 

o 

H 

1 

= c 

OO 

2 b. sin co. u + u 

J 3 


X 



1- J 


Lj=i J 


(9) 


The form of these equations suggests that by fitting the intensity at 2n 
channels we may obtain the unique set of n Fourier sine terms of amplitude 
bj and frequency prescribed by the observations. 

Let us check the method with a quadratic case, seeing to what extent 
we can reconstruct the familiar profile 


B(u) - 1 - exp (-u) 
on the basis of the four intensity observations 


l(0,lA i )-B(0) =a ± = f e" U dB(u) 


x i+ l 


1 , 2 , 3 , 4 . 


( 10 ) 


( 11 ) 
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Substitution of these values into Eq. (9) yields the following four 
relations which must be satisfied 


1 bl^l b 2 u 2 

2 = 2 + 2 

l+wq l+w 2 


1 

3 


2bi<i)i 

4+Mq 


+ 


2h 2 u 2 

*4 


1 

4 


^hj^l 3b 2 2 

2 + 2 
9+Wq 9+^ 


^l^l + llb 2 U} 2 

16+^ 2 l6+w 2 2 


( 12 ) 


After clearing of fractions and eliminating b and b 2 from the equations, 
we obtain the following two equations which must be satisfied 


with 


c - c 
1 2 

( 2+ 2 ) 2 2 
K w 2 j 3 1 w 2 

- 0 

°4 “ °5 

( 2+ 2 ) 4- 2 2 
H “2 6 W 1 w 2 

= 0 

c = 

27« q - 64 a_ + a = 

_1Z 

1 

3 2 1 

180 


5 5 3 


c = 
2 

3a^ + 
5 

l6a 

2 

15 

i 

^!h r 

ii 

_JZ. 

180 

c — 
3 

15 

4a 2 

15 

+ 

ii 

J=Z 

180 

c = 
4 

64q^ — 
7 

324 otj 
35 

+ 8a — 
2 

5 

1 

21 

c = 
5 

-4 a + 

k 

7 

36a 

35 

-2a = 
2 

5 

1 

105 


(13) 


(14) 
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°6 = !i " ^ + fT iji 
28 35 10 84 


The characteristic equation must be a binomial with roots at co 2 and 2 

1 2 

(x - Wi 2 ) (x -o> 2 2 ) = x 2 - ( co 1 2 t W2 2 ) x + o. ( 15 ) 

This equation is compatilbe with the Eqs. (13) if and only if the fol- 
lowing determinant vanishes 


c 

1 


c 

4 


c 

2 


c 

5 



( 16 ) 


yielding as the characteristic equation 


x 2 - c l c 6 ~ C 3 C 4 


— GqC 


3 5 


c l c 5 ~ c 2 C A = 
C 2 c 6 -c 3 c 5 


x 2 


ii X + 120 
11 55 


= 0. 


(17) 


This equation possesses the roots 

w-j 2 = 0.29350082 , co = ± 0.54175716 

(18) 

co 2 = 7.43377190 , w - + 2.72649443. 

We determine b_^and b^ now by back substitution of these roots into any 
pair of Eqs. (12) 

b l = 1.10607838, b^ = 0.11364969 ( 19 ) 
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Thus we have inferred from the four intensity measurements the following 
two Fourier sine terms 


B(u) 


b sin o) u + b sin u> u 
11 2 2 


= 1.10607838 sin 0.541757l6u 
+ 0.11364969 sin 2.72649443u. 


( 20 ) 



Figure 1. Nonlinear Fourier Inversion (inferred 
values versus actual profile.) 


Figure 1 displays the inferred compared to the actual profile. The 
agreement is highly gratifying. Particularly encouraging is the fact that 
the lower frequency component has some eleven times the weight of the higher 
term, an indication that the method converges rapidly. Although the inten- 
sities were chosen at integral values of the absorption coefficient, this need 
not have been the case, since the inversion algorithm is independent of 
channel choice. Separation is desirable, however, to avoid near singular 
matrices. 

Although the algebra becomes tedious as more Fourier terms are inferred 
it should be emphasized that once the channel sites are fixed, the routine 
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leading to the coefficient of the characteristic equation can be fixed once 
and for all. A program for direct time computer readout is certainly 
accessible. 

The algorithm will surely respond to noise in the data by having one or 
more of the roots of the characteristic equation becoming negative. In fact 
one can expect that the algorithm not only will discriminate against noise 
but infer which one (or more) of the channels is in error. This could be 
considered the ultimate in inversion, a method which would not only infer all 
the valid information implicit in the observations but also which channels are 
in error and by how much. 

3.2 Algorithm Response to Noise 

Let us examine now the reaction of the nonlinear Fourier inversion to 
data noise. As a test case we have retained the intensities in three of the 
directions associated with the sample profile [Eq. (10)], perturbing the 
third channel alone. The intensities now read 


~ 1/2 

a 2 = 1/3 

a 3 = 1/3 

a, - 1/5. 


( 21 ) 


The straightforward application of the algorithm as before yields the 
perturbed characteristic equation [cf. Eq.(l7)] 


* 2 + * 


1656 

361 


= 0 , 


( 22 ) 


with the roots, one positive and one negative, 


a^ 2 = + 0.55472626 
w 2 2 = - 8.26940770. 


(23) 


We determine as well by back substitution the weights 


bu = °* 78193 720 

(24) 

b^c^ = 0.02138690. 
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Our algorithm has therefore inferred from the perturbed intensity 
measurements , the following interpolation formula for the intensity 


2 

1(0,1/*) =2 * h.i^.i 

j=l 2 2 

* + cuj 

= 0.78193 720 x + 0.02138690 x , 
k 2 -K). 55472626 x 2 - 8 . 26940770 


( 25 ) 


and the following temperature profile 


2 

B(u) =2 b sin u> u 

j=l J j 

= 1.04986 sin ( 0 .7448u) + 0.007437i sin (2.8757iu). 


( 26 ) 


These results deserve the closest attention and interpretation. We see 
from the intensity formula [Eq. (25)] that the algorithm has inferred a pole 
at l/x = 1 /cop = 0.34774 to fit the erroneous measurement at l/x = 1/3. The 
contribution of the negative root at l/x = 1/3 is 


I (0,1/3) = 0.02138690/3 = .08781. 

2 l-(8. 26940770/9) (27) 


The actual displacement is (l/3)-(l/4) = (//l2) = .08333, which is within % 
of the true value. Moreover we see that the negative root term affects only 
slightly the three valid channels. This high fitting specificity follows 
from the concentration of the amplitude excursion near the pole. 

Turning now to the inferred temperature profile, [Eq. (26)] we see that 
there is only one valid term. The second term is deemed inadmissible because 
of the imaginary character of the amplitude and frequency. We note, however, 
that the contribution of the second term is down some two orders of magnitude 
from the first. The valid term that remains preserves with commendable 
fidelity the character of the actual profile. 

The inadmissible root response of the inversion algorithm to data error 
is reminiscent of the spline function inversion. Physically this is evidence 
that there does not exist any Fourier sine pair of arbitrary amplitudes and 
frequencies which will yield the observed intensities [Eq. (21)]. Hence we 
must ascribe the negative root to some effect, such as data error, extrinsic 
to the atmospheric emitting sources. 


463 


towards an optimal inversion method for remote atmospheric sensing 


These results underlie the optimism we feel that the ideal inversion 
procedure is within grasp and the expectation with which we view its appli- 
cation to real data. 
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